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We present maximally-fast numerical algorithms for conserved coarsening systems that are stable 
and accurate with a growing natural time-step At — AtV^ . For non-conserved systems, only effec- 
tively finite timesteps are accessible for similar unconditionally stable algorithms. We compare the 
scaling structure obtained from our maximally-fast conserved systems directly against the standard 
fixed-timestep Euler algorithm, and find that the error scales as \/A — so arbitrary accuracy can 
be achieved. 
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Phase-ordering kinetics studies the evolution of struc- 
ture after a quench from a disordered phase into an or- 
dered phase. The later stages of most phase-ordering pro- 
cesses in simple systems show universal scaling behavior 
described by a single growing length scale which increases 
as a power law in time, L{t) t", where < a < 1 0. 
For the scalar order-parameter systems considered in this 
paper, a = 1/2 and 1/3 for non-conserved and conserved 
dynamics, respectively 0] . While these growth exponents 
and their universality can be understood in terms of in- 
terfacial motion leading to domain coarsening 0, the 
time-independent scaled structure that result is less well 
understood. 

Computer simulation is an effective technique to sys- 
tematically study these non-linear non-equilibrium coars- 
ening systems. To maintain accuracy, the discretized dy- 
namics must move interfaces at most a small fraction of 
the interfacial width, ^, in a single timestep At. This 
determines a maximal or natural timestep of coarsening 
systems, M^at ~ ^/(dL/dt) ~ t^"", that grows in time. 
Unfortunately common time-discretizations are unstable 
for timesteps above a fixed threshold determined by the 
lattice spacing Ax Any such fixed timestep algo- 
rithm is increasingly inefficient at late times compared 
to the natural timestep. Various algorithms have been 
proposed to make simulations more efficient, including 
the cell-dynamical-scheme and Fourier spectral meth- 
ods [Bj. However, these approaches still require a fixed 
time-step for numerical stability. 

There is a newly developed class of unconditionally 
stable semi- implicit algorithms 0,13 that impose no sta- 
bility constraints on the timestep At, which is then de- 
termined by accuracy considerations. Since we generally 
expect larger At to lead to larger errors, there is a trade- 
off between speed and accuracy. This tradeoff is best 
resolved by picking growth rates for At that induce an 
error in the correlations that is approximately constant 
in magnitude throughout the scaling regime, where the 
magnitude can be chosen to be less than other system- 
atic sources of error such as initial transients or finite-size 
effects. While errors of single growing timesteps can be 



small , this begs the question of how much those single- 
step errors accumulate in correlations at late times after 
the quench. For example, we might expect that some 
types of single-step errors would be benign, given the ir- 
relevance of small amounts of random thermal noise to 
the scaled structure 0. Nevertheless, we would expect 
timesteps growing faster than At„at to lead to unaccept- 
able levels of error. 

In this paper, we compare the scaled correlations of 
unconditionally stable "Eyre" dynamics driven with a 
growing timestep At with the correlations evolved with 
an explicit Euler update. The latter, while slow, provides 
an arbitrarily accurate reference at late times. For con- 
served Cahn-Hilliard dynamics 9], we find that the maxi- 
mal difference of the scaled structure has a constant mag- 
nitude in the scaling regime when the Eyre algorithms are 
driven at the natural or maximal timestep 



At = AtV^, 



(1) 



for conserved dynamics. Here tg is a natural "structural 
time" determined by the decreasing system energy (see 
Eq. lO below). We find that this correlation error scales 
as ^/A for small A, and so can be made arbitrary small 
while retaining a maximally fast conserved coarsening al- 
gorithm which corresponds to moving interfaces a finite 
fraction of the interfacial width in every timestep. With 
a similar class of unconditionally stable algorithms for 
non-conserved dynamics 0, we find that only a fixed 
factor speedup is possible compared to the Euler algo- 
rithm as measured by the structural time ts- With a 
Fourier-space analysis of the dynamics, we explain how 
the natural timestep can accurately be used in conserved 
coarsening while only a fixed timestep is possible for non- 
conserved coarsening. 

Cahn-Hilliard dynamics of a conserved scalar field 
(j){r, t) are 

d(l)/dt^W'^6F/6(f>^-W^{(p + W^(p-<j>^), (2) 



where F = J d'^x[{V<j)y 



1)^/4] is the free en- 



ergy in d spatial dimensions with a double-well potential 
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corresponding to two distinct ordered phases at (j) = ±1. 
These dynamics can be semi-imphcitly discretized in time 
by 

^t+At + (1 - ai)AtV^4>t+At + (1 - a2)AtV'^i^t+At 

^(pt- MV\ai4>t + a2V^4>t - 0?), (3) 

where the discretized dynamics are unconditionally sta- 
ble for any At > when ai > 2 and 02 < 0.5 This 
equation implicitly defines the updated field, (j^t+At, and 
can be directly solved in Fourier space to give an Euler- 
like update 

Mt + ^t)^<Pk{t) + At,ff{k)^k (4) 

where the /c-dependent effective timestep is 

Ateff{k) = At/{1~ AtK), (5) 

where Ak is the Fourier-transformed Laplacian (Ak = 
—k^ in the continuum limit), and 

if = (ai - l)Ak + (02 - 1)A^ (6) 

We note that K < for unconditionally stable algo- 
rithms. 

Our numerical work is done in two-dimensions (2d), 
with systems of linear size Loo = 256 (at least 200 sam- 
ples) and Loo = 512 (at least 20 samples). We use a 
lattice spacing Ax — 1 and periodic boundary condi- 
tions. For Euler discretizations of the conserved dynam- 
ics Eq. Q with a\ = a2 = 1, we use At = 0.03, while for 
unconditionally stable discretizations we use ai = 3 and 
02 = 0. 

The fc-dependent effective timestep in Eq. Q led us 
to investigate the actual timestep taken for a given algo- 
rithmic timestep At. To do this we exploit the power-law 
decay of the free-energy density, e = F/V , in the late- 
time scaling regime, e ~ 1/L ^ t^" ^J, to introduce the 
"structural" time 

ts = i?e"i/", (7) 

where a = 1/3 or 1/2 for conserved or non-conserved 
coarsening, respectively, and B is chosen so that Ats — 
At for small At in the late-time scaling regime. For 
conserved dynamics we use B = 0.286, while for non- 
conserved Bnc = 0.105. The evolution of the structural 
time allows us to measure the real speed-up of our coars- 
ening algorithms, and it will also provide important in- 
sight into errors of the scaled structure. In general we 
find Ats/ At < 1 (see Fig. |21 and O below), and so we 
conservatively use the structural time ts to drive our al- 
gorithms, as in Eq. (Q). 

We find that using the natural timestep, Eq. Q), 
leads to accurate correlations in the scaling regime — as 
compared to systems evolved with a simple Euler time- 
discretization. We measure S{k,t) = {(j)k4'-k), where the 




FIG. 1: The average scaled structure S{x) = e^S{xe) vs. x = 
k/e for Loo ~ 512 system at ts — 1024 with the Euler update 
(circles) with At = 0.03 and the Eyre update ("-I-") with 
At — At^^^ where A — 0.01. Triangles indicate the absolute 
difference between Eyre and Euler updates, AS{x). 



angle-brackets indicate an average over orientations and 
initial conditions. We obtain the scaling form using the 
energy density, e, so that S{x) = e^S{xe) is a scaling 
function of a; = k/e. In Fig. ^ we plot S{x) vs. x to il- 
lustrate the excellent overlap between Euler (circles) and 
Eyre ("-F") dynamics with At = AtV^ and A = 0.01. 
To quantify the error we take the absolute value of the 
maximal value of the difference between the scaled struc- 
tures (shown for general k/e with triangles). We find 
that this maximum difference is approximately constant 
in magnitude throughout the scaling regime. 

We average this absolute error over the scaling regime, 
as determined by the scaling collapse of the scaled struc- 
ture. We observe a small A-dependent difference between 
the Euler algorithm and the naturally driven Eyre algo- 
rithm, as shown in Fig.|21 By repeating the measurement 
for different system sizes, both L^o = 256 (open circles, 
averaged over 4 times in tg G [60,190]) and Loo = 512 
(open squares, averaged over 9 times in ts £ [60, 1500]), 
we confirm that finite-size effects are not significant. We 
have enough independent Euler samples that no baseline 
errors due to residual stochastic effects are seen — the er- 
rors shown are the systematic error due to A. We observe 
an approximately \/A dependence on the average error. 
This implies that arbitrarily accurate measurements of 
the scaled structure can be made with maximally driven 
Eyre algorithms. 

From the effective timestep, Eq. ((SI, we can see quali- 
tatively why the natural timestep, Eq. is accurate for 
conserved updates. Since Fourier modes k 1/L ^ t^^^^ 
correspond to the domain structure, then for the natu- 
ral timestep with very small A the effective timestep is 
only /c-dependent deep in the Porod tail where kL ^ 1. 
However, the Porod tail jl| is simply the reflection of the 
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FIG. 2: The Eyre structural error vs. A for Eyre algorithms 
driven at the natural timestep At — At^^^ . The error is the 
maximum absolute difference between the Eyre correlations 
and the Euler correlations. Shown is data for Loo ~ 256 (open 
squares) and Loo ~ 512 (open circles). Also shown with filled 
symbols are the asymptotic values of 1 — Ats/At for various 
A. Both errors exhibit \/A behavior, as discussed in the text. 



amount of interface in the system at kL « 1. Systematic 
errors will arise as A gets larger and the /c-dependence of 
Ateff becomes more significant at kL « 1. 

Previous real-space single-step error analysis indi- 
cated that the natural timestep would not be accurate 
for direct Eyre algorithms like Eq. (PJ because of growing 
errors close to interfaces 8j. We have examined interfa- 
cial profiles of naturally driven algorithms, as probed by 
the k ~ tail of the unsealed structure, and found 
quantitative agreement with Euler algorithm (see large 
X data in Fig. Single-step errors must then corre- 
spond to errors in the interfacial motion, which can be 
probed by considering the difference between the algo- 
rithmic timestep At and the resulting structural timestep 
Atg. We find that Atg/At < 1 in general, which indi- 
cates reduced interfacial motion. In Fig. [21 we show with 
filled symbols that 1 — Ats/At scales as ^/A for small A. 
We recover this asymptotic result using an energy-scaling 
argument 1, 2] in general spatial dimension d. 

From Eq. 0, with a = 1/3 for conserved dynamics, 
we have e — B^^^t^^^^ so that 

Ats = -iB-^/'^Aetl'^. (8) 

On the other hand we can integrate the energvdissipated 
in one timestep for each Fourier component [2| , 

Ae ~ I d''k/{2TTf{{SF/Scf>k)Aq^k), (9) 

= - j d''k/{2^fk-^At,ff{k,At)Tk (10) 

where the second line uses the time-derivative 4>-k = 
-k^SF/S(t>k from Eq. 0, A0fc = Ateff(j>k from Eq. igj). 



and the time-derivative scaling function 

Tk = i'i^k^-k) = L'L^-^hikL) (11) 

where the scaling form is shown. Tk is expected to have 
a Porod-like h{x) ~ x^~'^ tail for x ^ 1 2\ and this has 
been observed in d = 2 ^3 ■ With these asymptotics and 
Ateff from Eq. ((SJ , the Ae integral converges and hence 
becomes time-independent as the UV cutoff 0{L/^) be- 
comes large. For At = Ats we obtain 

Ats/At (X / dxx'^-^h{x)/{l + A'x"^) (12) 
Jo 

where A' = A{ai — 1)/Lq and L = Lotl^"^ in the scaling 
regime. For small A, the leading contribution is 0{VA') 
from the large x regime. Since Atg/At = 1 in the limit 
of Ai ^ when ^ = 0"*" , we have 

1 - Ats/At cx VA' + 0{A') (13) 

This behavior is observed, as shown by the filled points 
in Fig. El 

What is the connection, if any, between the timestep 
error 1 — Ats/At and the structural error, both of which 
exhibit ^/A dependence at small A for natural timesteps? 
As discussed before, 1 — Atg/At ~ ^/A > indicates an 
error of (reduced) interfacial motion. It is reasonable that 
this error shows up in the correlations at the same order, 
0{^/A). It is interesting that this error accumulates into 
a constant contribution to the scaled correlations within 
the scaling regime. 

We now consider non-conserved AUen-Cahn coarsening 
dynamics, which are governed by 

= -(5F/(50 = + V2(/)-03. (14) 

These dynamics can be semi-implicitly discretized in time 
by 

4>t+At + (ai - l)At4>t+At + (a2 - l)AtW^4>t+At 

= (j)t + At{ai(t>t+a2V^<j)t- (/yt), (15) 

where the discretized dynamics are unconditionally sta- 
ble for any At > when ai > 2 and 0-2 < 0.5 7J. In the 
same spirit as conserved dynamics, we obtain an effective 
timestep 

Ateffik,At) ^ At/{1~ AtN), (16) 

where N = (1 — ai) + (1 — 02)Ak and < for stable 
algorithms with ai > 2 and 02 < 0.5 01 . This directly 
implies that At^ff < 1/(1 — ai) even when At — > 00, 
so that this class of unconditionally stable non-conserved 
algorithms effectively cannot be accelerated. We confirm 
this by calculating and measuring Ats when At = 00. 

We adapt the development in Eq. JSJ - H12(l for non- 
conserved dynamics. We have e = BlJ^ts so that 

Ats = -2B-y'Aetl/\ (17) 
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FIG. 3: For non-conserved update driven with At = 
oo, we plot Ats = Atsy^ (tti — 1)(1 — 02) vs. a = 

— 02)/ (ffli — 1) showing the prediction (solid hne with 
5 = 0.85)) of Eq. CT . 



Integrating the energy dissipated in one timestep, 

Ae~- J d'^k{2TTy'^Ateff{k, At)Tk (18) 

where we use (f)-k — ~SF/S(t)k- Using Eq. (|16|l and 
Eq. 1)11(1 with L = L^t^J^ and Aj. = — fc^ we can solve 
for At, when At = 00: 



At, 




x"^ ^h{x)dx 



(ai - l)L2 + a;2(l-a2) 



(19) 



l)(l-a2)' 



(20) 



where we take the late-time L —* 00 limit in the second 
line and have used hnc{x) ~ x^"'' for x 3> 1 [ll 1^ Hof. 
without which Atg is not time-independent. The overall ^ 
factor on the second line comes from imposing Atg — At 
for small At to determine Bncj before the At — > 00 
limit is taken. We find a constant Ats that depends 
only on ai and 02, as well as the inverse UV cutoff ^. 
In_ Fig. for 0? = 2, we plot the measured asymptotic 
Ats = Ats-\/(ai — 1)(1 — 02) (averaged over the time- 
independent scaling regime with variances shown) vs. 
a = a/(1 — a2)/{ai — 1). We observe data collapse for 
a variety of oi and 02, as suggested by Eq. (|20|l . How- 
ever we only agree with the calculated functional form 
for small a - the solid line indicates best fit by eye 
with ^ = 0.85. The discrepancy appears to be due to 
a dependence of the time-derivative correlations (data 
not shown) , indicating significant systematic errors when 
At — !■ 00 despite the finite effective timestep. 

In summary, we have shown that for conserved dynam- 
ics we obtain accurate scaled correlations for maximally 
fast algorithms with At = At'l^^ . The structural error 



behaves as vG4 for small A. These results are consis- 
tent with previous real-space single-step error-analysis 
away from interfacial regions. For these maximally fast 
algorithms the relative speedup with respect to a fixed 
timestep is of the order (Lac/Ax)'^, where Loo /Ax is the 
discretized linear system size. Maximally fast algorithms 
provide the most efficient means to reach the scaling limit 
for large systems. 

Surprisingly, a similar class of algorithms does not lead 
to acceleration of non-conserved dynamics. Only a con- 
stant time-step is observed, as measured by the structural 
time, even when the unconditionally stable algorithm is 
driven with At ~^ 00. 

We expect that these unconditionally-stable semi- 
implicit algorithms find broader application. In the sys- 
tems investigated so far, the regime of unconditional 
stability coincides with the easily determined regime of 
linear-stability. For coarsening systems that exhibit a 
growing natural timestep, large accelerations are possi- 
ble with growing At. While the non-conserved dynam- 
ics cautions us that effective acceleration is not always 
produced, it also illustrates the diagnostic value of the 
ratio Ats/ At. Indeed, 1 — Ats/ At appears to provide a 
good proxy for systematic structural errors in the scaling 
regime. 
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